Generalized linear mixed effects models were appropriate for this analysis because the data contain sources of non-independence. Individual plants were located in the same site as other individuals and were likely genetically similar. Therefore, we considered site a random effect.
The response variable is total seeds, which is discrete and non-zero. The interaction between deviation and the early/late factor was a fixed effect, as well as the interaction between plot treatment and conspecific density. We would have included soil moisture as a covariate and fixed effect in our model, but soil moisture was not important for seed production. Each species is analyzed separately.
Mertensia
We started by limiting the individuals to the open treatment (excluding the hand-pollinated treatment) and by checking the distribution of the total seed counts.
mert.open<-mert[(mert$treat=="open"),] #removed hand pollination treatment
mert.nbinom<-fitdist(mert.open$dev.seed, "nbinom") #fitting Mertensia data to negative binomial distribution using fitdistplus
plot(mert.nbinom) #plotting to see the fit
The Mertensia seed data fits a negative binomial distribution. The glmmTMB package is a good package for modeling data with a negative binomial. The glmmTMB package is also suitable for data that could be overdispersed or zero-inflated.
Main Model
#mert.glmmtmb<-glmmTMB(dev.seed ~ flowers:(deviation*early.late) + flowers:number.conspecifics/plot.treat + (1|site), family = "nbinom2", data = mert.open) #putting count data into glmmTMB for DHARMa package
#summary(mert.glmmtmb) #summary of glmmTMB model
#mert.counts.simulation<-simulateResiduals(fittedModel = mert.glmmtmb) #creating simulation for Mertensia mixed effects model with DHARMa package
#plot(mert.counts.simulation) #plotting simulation
#testDispersion(mert.counts.simulation) #checking for overdispersion
#testZeroInflation(mert.counts.simulation) #checking for zero inflation
We detected a significant effect of early vs. late blooming on the total number of seeds in Mertensia individuals, but when we checked for overdispersion and zero-inflation in the data using the DHARMa package, the data were slightly zero-inflated. To address this, we adjusted our glmmtmb model to include zero-inflation.
mert.glmmtmb<-glmmTMB(dev.seed ~ flowers:(deviation*early.late) + flowers:number.conspecifics/plot.treat + (1|site), ziformula=~1, family = "nbinom2", data = mert.open) #putting count data into glmmTMB for DHARMa package
summary(mert.glmmtmb) #summary of glmmTMB model
## Family: nbinom2 ( log )
## Formula:
## dev.seed ~ flowers:(deviation * early.late) + flowers:number.conspecifics/plot.treat +
## (1 | site)
## Zero inflation: ~1
## Data: mert.open
##
## AIC BIC logLik deviance df.resid
## 386.2 402.4 -184.1 368.2 36
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.07728 0.278
## Number of obs: 45, groups: site, 4
##
## Overdispersion parameter for nbinom2 family (): 3.03
##
## Conditional model:
## Estimate Std. Error z value
## (Intercept) 2.613e+00 3.761e-01 6.949
## flowers:deviation 2.433e-03 1.540e-03 1.579
## flowers:early.latelate 1.931e-02 1.009e-02 1.914
## flowers:number.conspecifics -2.368e-05 8.772e-05 -0.270
## flowers:deviation:early.latelate -2.693e-03 2.934e-03 -0.918
## flowers:number.conspecifics:plot.treatmanip 2.973e-04 1.937e-04 1.535
## Pr(>|z|)
## (Intercept) 3.68e-12 ***
## flowers:deviation 0.1143
## flowers:early.latelate 0.0557 .
## flowers:number.conspecifics 0.7872
## flowers:deviation:early.latelate 0.3587
## flowers:number.conspecifics:plot.treatmanip 0.1248
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.3800 0.5513 -4.317 1.58e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mert.counts.simulation<-simulateResiduals(fittedModel = mert.glmmtmb) #creating simulation for Mertensia mixed effects model with DHARMa package
plot(mert.counts.simulation) #plotting simulation
testDispersion(mert.counts.simulation) #checking for overdispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.90351, p-value = 0.784
## alternative hypothesis: two.sided
testZeroInflation(mert.counts.simulation) #checking for zero inflation
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0695, p-value = 1
## alternative hypothesis: two.sided
When number of flowers is taken into account, seed set for Mertensia individuals may be affected by blooming before or after the population peak. Zero-inflation is accounted for in the new model, and the model for Mertensia count data is no longer zero-inflated.
To figure out the direction of the effect, we looked at the plots (shown below). Blooming after the population peak may increase seed set for Mertensia individuals.
Pollen Limitation Model
To test the pollen limitation individuals, we created a separate fixed effects model. This model has site and plot treatment as random effects and the interaction between pollen limitation treatment and number of flowers as a fixed effect. We included this interaction because the difference in seed set between the open and hand treatments is heavily dependent on the number of flowers. The interaction accounts for the wide variation of number of flowers and ultimately developed seeds. The response variable is still total number of developed seeds in Mertensia individuals.
#mert.treat.glmmtmb<-glmmTMB(dev.seed ~ treat/flowers + (1|site/plot.treat), family = "nbinom2", data = mert) #Mertensia model for negative binomial data with glmmTMB package
#summary(mert.treat.glmmtmb) #summary of model output
#mert.counts.treat.simulation<-simulateResiduals(fittedModel = mert.treat.glmmtmb) #creating simulation for Mertensia mixed effects model with DHARMa package
#plot(mert.counts.treat.simulation) #plotting simulation
#testDispersion(mert.counts.treat.simulation) #checking for overdispersion
#testZeroInflation(mert.counts.treat.simulation) #checking for zero inflation
We detected an effect of pollen limitation treatment on Mertensia total developed seeds, but we had to test the model assumptions of our pollen limitation model as well. The data in the pollen limitation model for Mertensia are not overdispersed, but they are zero-inflated. We then adjusted our model to include zero-inflation and ran the overdispersion and zero-inflation tests again.
mert.treat.glmmtmb<-glmmTMB(dev.seed ~ treat/flowers + (1|site/plot.treat), ziformula=~1, family = "nbinom2", data = mert) #Mertensia model for negative binomial data with glmmTMB package
summary(mert.treat.glmmtmb) #summary of model output
## Family: nbinom2 ( log )
## Formula: dev.seed ~ treat/flowers + (1 | site/plot.treat)
## Zero inflation: ~1
## Data: mert
##
## AIC BIC logLik deviance df.resid
## 783.1 803.7 -383.6 767.1 89
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## plot.treat:site (Intercept) 0.02665 0.1633
## site (Intercept) 0.07150 0.2674
## Number of obs: 97, groups: plot.treat:site, 8; site, 4
##
## Overdispersion parameter for nbinom2 family (): 5.6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.376715 0.242589 9.797 < 2e-16 ***
## treatopen-hand -0.097678 0.255614 -0.382 0.702
## treatopen:flowers 0.024074 0.004401 5.470 4.51e-08 ***
## treatopen-hand:flowers 0.026340 0.004167 6.320 2.61e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.9293 0.4668 -6.275 3.49e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mert.counts.treat.simulation<-simulateResiduals(fittedModel = mert.treat.glmmtmb) #creating simulation for Mertensia mixed effects model with DHARMa package
plot(mert.counts.treat.simulation) #plotting simulation
testDispersion(mert.counts.treat.simulation) #checking for overdispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.85214, p-value = 0.592
## alternative hypothesis: two.sided
testZeroInflation(mert.counts.treat.simulation) #checking for zero inflation
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0121, p-value = 1
## alternative hypothesis: two.sided
The data are no longer zero-inflated, and we still detected an effect of pollen limitation on the Mertensia seed counts. Mertensia individuals are pollen limited.
Figures
Below is the plot of seed set vs. relative position of blooming in the population for Mertensia individuals. Plot treatment is indicated by individual color.
ggplot(mert.open, aes(x=mert.open$relative.position, y=mert.open$dev.seed, group=mert.open$plot.treat)) +
labs(y="Total Developed Seeds",x="Days Relative to Population Peak Bloom", title="Effect of Mertensia Blooming Time on Total Developed Seeds") +
geom_point(aes(color=plot.treat),size=3, alpha = 0.8) +
scale_color_brewer(palette="Dark2") + labs(color="Plot Treatment") +
scale_color_manual(labels= c("control","accelerated"), values = c("control"="darkred","manip"="darkorange")) +
geom_smooth()
ggplot(mert.open, aes(x=mert.open$relative.position, y=mert.open$dev.seed, group=mert.open$early.late)) +
labs(y="Total Developed Seeds",x="Days Relative to Population Peak Bloom", title="Effect of Mertensia Blooming Time on Total Developed Seeds") +
theme_classic() +
geom_point(aes(color=early.late),size=3, alpha = 0.8) +
scale_color_brewer(palette="Dark2") + labs(color="Before or After Peak") +
scale_color_manual(labels= c("before peak","after peak"), values = c("early"="darkred","late"="darkorange")) +
geom_smooth(method = "glm", aes(color=early.late))
We then plotted the interaction between deviation and early/late instead of relative position in order to fit a linear model.
ggplot(mert.open, aes(x=mert.open$deviation, y=mert.open$dev.seed/mert.open$flowers, group= mert.open$early.late)) +
labs(y="Total Developed Seeds per Flower",x="Days Relative to Population Peak Bloom", title="Effect of Mertensia Deviation from Peak on Total Developed Seeds") + theme_classic() +
theme_classic() +
geom_point(size=3, alpha = 0.8, aes(color=early.late)) +
geom_smooth(method="glm",aes(color=early.late)) +
scale_color_brewer(palette="Dark2") +
labs(color="Before or After Peak") +
scale_color_manual(labels= c("before peak","after peak"), values = c("early"="darkred","late"="darkorange"))
Then we plotted the effect of conspecific density and plot treatment on total developed seeds.
ggplot(mert.open, aes(x=mert.open$number.conspecifics, y=mert.open$dev.seed, group= mert.open$plot.treat)) +
labs(y="Total Developed Seeds",x="Conspecific Density", title="Effect of Mertensia Conspecific Density on Total Developed Seeds") +
theme_classic() +
geom_point(size=3, alpha = 0.8, aes(color=plot.treat)) +
geom_smooth(method="glm",aes(color=plot.treat)) +
scale_color_brewer(palette="Dark2") +
labs(color="Plot Treatment") +
scale_color_manual(labels= c("control","accelerated snowmelt"), values = c("control"="darkgreen","manip"="blue"))
Finally, we plotted all of the data with the pollinator active period. The pollinator active period was the time during which we saw interactions between pollinators and our species. The active period is shaded.
mert$midpoint2<-format(as.Date(mert$midpoint), format="%m-%d")
ggplot(mert, aes(x=mert$midpoint2, y=mert$dev.seed, group= mert$treat)) +
labs(y="Total Developed Seeds",x="Date", title="Effect of Mertensia Blooming Time and Pollinator Activity on Total Developed Seeds") +
theme_classic() +
geom_point(size=3, alpha = 0.8, aes(color=treat)) +
geom_smooth(aes(color=treat)) +
scale_color_brewer(palette="Dark2") +
labs(color="Open or Hand-pollinated") +
scale_color_manual(labels= c("open","open-hand"), values = c("open"="darkred","open-hand"="darkorange")) +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
annotate("rect",xmin="06-19", xmax="07-03",ymin=0,ymax=Inf,alpha=0.1)
Delphinium
For Delphinium, we again limited the individuals to open treatment individuals.
delph.open<-delph[(delph$treat=="open"),] #removed hand pollination treatment
delph.nbinom<-fitdist(delph.open$total.seeds, "nbinom") #fitting negative binomial distribution using fitdistplus package
plot(delph.nbinom) #plotting distribution fit
The Delphinium data fits a negative binomial distribution. We used the glmmTMB package again.
Main Model
delph.glmmtmb<-glmmTMB(total.seeds ~ X.flowers:(deviation*early.late) + X.flowers:number.conspecifics/plot.treat + (1|site), family = "nbinom2", data = delph.open) #Delphinium model for negative binomial data with glmmTMB package
#removed treat fixed effect (only open)
summary(delph.glmmtmb) #summary of model output
## Family: nbinom2 ( log )
## Formula:
## total.seeds ~ X.flowers:(deviation * early.late) + X.flowers:number.conspecifics/plot.treat +
## (1 | site)
## Data: delph.open
##
## AIC BIC logLik deviance df.resid
## 403.9 417.6 -193.9 387.9 33
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.02556 0.1599
## Number of obs: 41, groups: site, 4
##
## Overdispersion parameter for nbinom2 family (): 1.37
##
## Conditional model:
## Estimate Std. Error z value
## (Intercept) 2.990040 0.369707 8.088
## X.flowers:deviation 0.026479 0.011771 2.249
## X.flowers:early.latelate 0.137212 0.168217 0.816
## X.flowers:number.conspecifics 0.001297 0.001946 0.666
## X.flowers:deviation:early.latelate -0.010662 0.036032 -0.296
## X.flowers:number.conspecifics:plot.treatmanip -0.002071 0.001801 -1.150
## Pr(>|z|)
## (Intercept) 6.09e-16 ***
## X.flowers:deviation 0.0245 *
## X.flowers:early.latelate 0.4147
## X.flowers:number.conspecifics 0.5051
## X.flowers:deviation:early.latelate 0.7673
## X.flowers:number.conspecifics:plot.treatmanip 0.2501
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
delph.counts.simulation<-simulateResiduals(fittedModel = delph.glmmtmb) #creating simulation for Delphinium mixed effects model with DHARMa package
plot(delph.counts.simulation) #plotting simulation
testDispersion(delph.counts.simulation) #checking for overdispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.81708, p-value = 0.496
## alternative hypothesis: two.sided
testZeroInflation(delph.counts.simulation) #checking for zero inflation
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 4.5872, p-value = 0.152
## alternative hypothesis: two.sided
When flower number is taken into account, deviation from the population peak has a significant effect on Delphinium developed seed count. We tested the model assumptions as we did with Mertensia seed count data. The model for the Delphinium count data is not overdispersed or zero-inflated.
To determine the direction of deviation effect, we looked at the plot. When flower number is taken into account, seed set of Delphinium individuals increase with deviation from the population peak.
Pollen Limitation Model
Again, we tested the pollen limitation in a separate mixed effects model.
#delph.treat.glmmtmb<-glmmTMB(total.seeds ~ treat/X.flowers + (1|site/plot.treat), family = "nbinom2", data = delph) #Mertensia model for negative binomial data with glmmTMB package
#summary(delph.treat.glmmtmb) #summary of model output
#delph.counts.treat.simulation<-simulateResiduals(fittedModel = delph.treat.glmmtmb) #creating simulation for Delphinium mixed effects model with DHARMa package
#plot(delph.counts.treat.simulation) #plotting simulation
#testDispersion(delph.counts.treat.simulation) #checking for overdispersion
#testZeroInflation(delph.counts.treat.simulation) #checking for zero inflation
We detected an effect of pollen limitation on the total developed seeds of Delphinium individuals. We again checked for overdispersion and zero-inflation in the pollen limitation data. The Delphinum pollen limitation model is overdispersed and zero-inflated. We adjusted our model to include the zero-inflation formula, and re-tested for overdispersion and zero-inflation.
delph.treat.glmmtmb<-glmmTMB(total.seeds ~ treat/X.flowers + (1|site/plot.treat), ziformula=~1, family = "nbinom2", data = delph) #Mertensia model for negative binomial data with glmmTMB package
summary(delph.treat.glmmtmb) #summary of model output
## Family: nbinom2 ( log )
## Formula: total.seeds ~ treat/X.flowers + (1 | site/plot.treat)
## Zero inflation: ~1
## Data: delph
##
## AIC BIC logLik deviance df.resid
## 937.0 957.5 -460.5 921.0 88
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## plot.treat:site (Intercept) 0.003478 0.05897
## site (Intercept) 0.071347 0.26711
## Number of obs: 96, groups: plot.treat:site, 7; site, 4
##
## Overdispersion parameter for nbinom2 family (): 1.97
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.54661 0.36803 6.920 4.53e-12 ***
## treatopen-hand 0.55500 0.47352 1.172 0.24117
## treatopen:X.flowers 0.29208 0.06842 4.269 1.97e-05 ***
## treatopen-hand:X.flowers 0.18467 0.06693 2.759 0.00579 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.2789 0.5892 -5.565 2.63e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
delph.counts.treat.simulation<-simulateResiduals(fittedModel = delph.treat.glmmtmb) #creating simulation for Delphinium mixed effects model with DHARMa package
plot(delph.counts.treat.simulation) #plotting simulation
testDispersion(delph.counts.treat.simulation) #checking for overdispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.73169, p-value = 0.096
## alternative hypothesis: two.sided
testZeroInflation(delph.counts.treat.simulation) #checking for zero inflation
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0604, p-value = 1
## alternative hypothesis: two.sided
The new model is not overdispersed or zero-inflated. We still see an effect of pollen limitation on the Delphinium count data. Delphinium individuals are pollen limited.
Figures
Below are the plots for the Delphinium individuals.
ggplot(delph.open, aes(x=delph.open$relative.position, y=delph.open$total.seeds, group=plot.treat)) +
labs(title="Effect of Delphinium Blooming Time on Total Developed Seeds", y="Total Developed Seeds",x="Days Relative to Population Peak Bloom") +
geom_point(aes(color=plot.treat),size=3, alpha = 0.8) +
scale_color_brewer(palette="Dark2") + labs(color="Plot Treatment") +
scale_color_manual(labels= c("control","accelerated"), values = c("control"="darkred","manip"="darkorange")) +
geom_smooth()
ggplot(delph.open, aes(x=delph.open$relative.position, y=delph.open$total.seeds, group=delph.open$early.late)) +
labs(y="Total Developed Seeds",x="Days Relative to Population Peak Bloom", title="Effect of Delphinium Blooming Time on Total Developed Seeds") +
theme_classic() +
geom_point(aes(color=early.late),size=3, alpha = 0.8) +
scale_color_brewer(palette="Dark2") + labs(color="Before or After Peak") +
scale_color_manual(labels= c("before peak","after peak"), values = c("early"="darkred","late"="darkorange")) +
geom_smooth(method = "glm", aes(color=early.late))
ggplot(delph.open, aes(x=delph.open$deviation, y=(delph.open$total.seeds/delph.open$X.flowers), group= delph.open$early.late)) +
labs(y="Total Developed Seeds per Flower",x="Days Relative to Population Peak Bloom", title="Effect of Delphinium Deviation from Peak on Total Developed Seeds") +
theme_classic() +
geom_point(size=3, alpha = 0.8, aes(color=early.late)) +
geom_smooth(method="glm",aes(color=early.late)) +
scale_color_brewer(palette="Dark2") +
labs(color="Before or After Peak") +
scale_color_manual(labels= c("before peak","after peak"), values = c("early"="darkred","late"="darkorange"))
ggplot(delph.open, aes(x=delph.open$number.conspecifics, y=delph.open$total.seeds, group= delph.open$plot.treat)) +
labs(y="Total Developed Seeds",x="Conspecific Density", title="Effect of Delphinium Conspecific Density on Total Developed Seeds") +
theme_classic() +
geom_point(size=3, alpha = 0.8, aes(color=plot.treat)) +
geom_smooth(method="glm",aes(color=plot.treat)) +
scale_color_brewer(palette="Dark2") +
labs(color="Plot Treatment") +
scale_color_manual(labels= c("control","accelerated snowmelt"), values = c("control"="darkgreen","manip"="blue"))
Finally, we plotted the pollinator active period with the total number of developed seeds and midpoint of individuals.
delph$midpoint2<-format(as.Date(delph$midpoint), format="%m-%d")
ggplot(delph, aes(x=delph$midpoint2, y=delph$total.seeds, group= delph$treat)) +
labs(y="Total Developed Seeds",x="Date", title="Effect of Delphinium Blooming Time and Pollinator Activity on Total Developed Seeds") +
theme_classic() +
geom_point(size=3, alpha = 0.8, aes(color=treat)) +
geom_smooth(aes(color=treat)) +
scale_color_brewer(palette="Dark2") +
labs(color="Open or Hand-pollinated") +
scale_color_manual(labels= c("open","open-hand"), values = c("open"="darkred","open-hand"="darkorange")) +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
annotate("rect",xmin="06-27", xmax="07-18",ymin=0,ymax=Inf,alpha=0.1)
Potentilla
We began by limiting the data to the open treatment Potentilla individuals.
pot.open<-pot[(pot$treat=="open"),] #removed hand pollination treatment
pot.nbinom<-fitdist(pot.open$total.seeds, "nbinom") #fitting negative binomial distribution using fitdistplus package
plot(pot.nbinom) #plotting distribution fit
The Potentilla seed data fits a negative binomial distribution. We used the glmmTMB package again.
Main Model
#pot.glmmtmb<-glmmTMB(total.seeds ~ X.flowers:(deviation * early.late) + X.flowers:number.conspecifics/plot.treat + (1|site), family = "nbinom2", data = pot.open) #Potentilla model for negative binomial data with glmmTMB package
#summary(pot.glmmtmb) #summary of model output
#pot.counts.simulation<-simulateResiduals(fittedModel = pot.glmmtmb) #creating simulation for Potentilla mixed effects model with DHARMa package
#plot(pot.counts.simulation) #plotting simulation
#testDispersion(pot.counts.simulation) #checking for overdispersion
#testZeroInflation(pot.counts.simulation) #checking for zero inflation
The Potentilla data are not overdispersed, but they are zero-inflated. To fix this, we had to include zero inflation in our model and again check for zero-inflation.
pot.glmmtmb<-glmmTMB(total.seeds ~ X.flowers:(deviation * early.late) + X.flowers:number.conspecifics/plot.treat + (1|site), ziformula=~1, family = "nbinom2", data = pot.open) #new Potentilla model for negative binomial data with glmmTMB package
summary(pot.glmmtmb) #summary of model output
## Family: nbinom2 ( log )
## Formula:
## total.seeds ~ X.flowers:(deviation * early.late) + X.flowers:number.conspecifics/plot.treat +
## (1 | site)
## Zero inflation: ~1
## Data: pot.open
##
## AIC BIC logLik deviance df.resid
## 280.7 291.4 -131.4 262.7 15
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 3.537e-10 1.881e-05
## Number of obs: 24, groups: site, 4
##
## Overdispersion parameter for nbinom2 family (): 2.95
##
## Conditional model:
## Estimate Std. Error z value
## (Intercept) 3.657e+00 4.095e-01 8.930
## X.flowers:deviation 1.220e-02 4.421e-03 2.760
## X.flowers:early.latelate 7.605e-02 2.733e-02 2.783
## X.flowers:number.conspecifics 1.631e-05 3.236e-04 0.050
## X.flowers:deviation:early.latelate -1.423e-02 4.629e-03 -3.074
## X.flowers:number.conspecifics:plot.treatmanip 4.595e-04 3.227e-04 1.424
## Pr(>|z|)
## (Intercept) < 2e-16 ***
## X.flowers:deviation 0.00578 **
## X.flowers:early.latelate 0.00539 **
## X.flowers:number.conspecifics 0.95981
## X.flowers:deviation:early.latelate 0.00211 **
## X.flowers:number.conspecifics:plot.treatmanip 0.15444
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.138 1.024 -3.065 0.00218 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pot.counts.simulation<-simulateResiduals(fittedModel = pot.glmmtmb) #creating simulation for Potentilla mixed effects model with DHARMa package
plot(pot.counts.simulation) #plotting simulation
testDispersion(pot.counts.simulation) #checking for overdispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.8079, p-value = 0.52
## alternative hypothesis: two.sided
testZeroInflation(pot.counts.simulation) #checking for zero inflation
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0638, p-value = 1
## alternative hypothesis: two.sided
When we incorporated zero-inflation into the model, we found a significant effect of deviation from the population peak, blooming before vs after the peak, and the interaction between the two. The direction of the effect, however, must be determined from the plots (below).
After looking at the plots, we see that when we account for the number of flowers, deviation has a greater negative effect on the late blooming Potentilla individuals.
Pollen Limitation Model
Below is the mixed effects model for the pollen limitation treatments.
#pot.treat.glmmtmb<-glmmTMB(total.seeds ~ treat/X.flowers + (1|site/plot.treat), family = "nbinom2", data = pot) #Potentilla model for negative binomial data with glmmadmb package
#summary(pot.treat.glmmtmb) #summary of model output
#pot.counts.treat.simulation<-simulateResiduals(fittedModel = pot.treat.glmmtmb) #creating simulation for Potentilla mixed effects model with DHARMa package
#plot(pot.counts.treat.simulation) #plotting simulation
#testDispersion(pot.counts.treat.simulation) #checking for overdispersion
#testZeroInflation(pot.counts.treat.simulation) #checking for zero inflation
We detected a significant effect of the pollen limitation treatment on the total number of developed seeds for Potentilla individuals, but the Potentilla pollen limitation data were overdispersed and zero-inflated. We addressed the zero-inflation issue by again including the zero-inflation formula into our models.
pot.treat.glmmtmb<-glmmTMB(total.seeds ~ treat/X.flowers + (1|site/plot.treat), ziformula=~1, family="nbinom2",data = pot) #Potentilla model for negative binomial data with glmmadmb package
summary(pot.treat.glmmtmb) #summary of model output
## Family: nbinom2 ( log )
## Formula: total.seeds ~ treat/X.flowers + (1 | site/plot.treat)
## Zero inflation: ~1
## Data: pot
##
## AIC BIC logLik deviance df.resid
## 502.1 516.4 -243.0 486.1 36
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## plot.treat:site (Intercept) 6.676e-09 8.171e-05
## site (Intercept) 1.193e-02 1.092e-01
## Number of obs: 44, groups: plot.treat:site, 6; site, 4
##
## Overdispersion parameter for nbinom2 family (): 1.86
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.87623 0.56142 6.904 5.04e-12 ***
## treatopen-hand 0.61864 0.68954 0.897 0.3696
## treatopen:X.flowers 0.05534 0.03138 1.764 0.0778 .
## treatopen-hand:X.flowers 0.01014 0.02816 0.360 0.7189
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.6271 0.6051 -4.342 1.41e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pot.counts.treat.simulation<-simulateResiduals(fittedModel = pot.treat.glmmtmb) #creating simulation for Potentilla mixed effects model with DHARMa package
plot(pot.counts.treat.simulation) #plotting simulation
testDispersion(pot.counts.treat.simulation) #checking for overdispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.80175, p-value = 0.336
## alternative hypothesis: two.sided
testZeroInflation(pot.counts.treat.simulation) #checking for zero inflation
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0204, p-value = 1
## alternative hypothesis: two.sided
The data are no longer zero-inflated or overdispersed, and we still detected a marginally significant effect of pollen limitation in the Potentilla individuals. Potentilla individuals are slightly pollen limited.
Figures
Below are the plots for Potentilla individuals.
ggplot(pot.open, aes(x=pot.open$relative.position, y=pot.open$total.seeds, group=plot.treat)) +
labs(title="Effect of Potentilla Blooming Time on Total Developed Seeds", y="Total Developed Seeds",x="Days Relative to Population Peak Bloom") +
geom_point(aes(color=plot.treat),size=3, alpha = 0.8) +
scale_color_brewer(palette="Dark2") + labs(color="Plot Treatment") +
scale_color_manual(labels= c("control","accelerated"), values = c("control"="darkred","manip"="darkorange")) +
geom_smooth()
ggplot(pot.open, aes(x=pot.open$relative.position, y=pot.open$total.seeds, group=pot.open$early.late)) +
labs(y="Total Developed Seeds",x="Days Relative to Population Peak Bloom", title="Effect of Potentilla Blooming Time on Total Developed Seeds") +
theme_classic() +
geom_point(aes(color=early.late),size=3, alpha = 0.8) +
scale_color_brewer(palette="Dark2") + labs(color="Before or After Peak") +
scale_color_manual(labels= c("before peak","after peak"), values = c("early"="darkred","late"="darkorange")) +
geom_smooth(method = "glm", aes(color=early.late))
ggplot(pot.open, aes(x=pot.open$deviation, y=pot.open$total.seeds/pot.open$X.flowers, group= pot.open$early.late)) +
labs(y="Total Developed Seeds per Flower",x="Days Relative to Population Peak Bloom",title="Effect of Potentilla Deviation from Peak on Total Developed Seeds") +
theme_classic() +
geom_point(size=3, alpha = 0.8, aes(color=early.late)) +
geom_smooth(method="glm",aes(color=early.late)) +
scale_color_brewer(palette="Dark2") +
labs(color="Before or After Peak") +
scale_color_manual(labels= c("before peak","after peak"), values = c("early"="darkred","late"="darkorange"))
ggplot(pot.open, aes(x=pot.open$number.conspecifics, y=pot.open$total.seeds, group= pot.open$plot.treat)) +
labs(y="Total Developed Seeds",x="Conspecific Density", title="Effect of Potentilla Conspecific Density on Total Developed Seeds") +
theme_classic() +
geom_point(size=3, alpha = 0.8, aes(color=plot.treat)) +
geom_smooth(method="glm",aes(color=plot.treat)) +
scale_color_brewer(palette="Dark2") +
labs(color="Plot Treatment") +
scale_color_manual(labels= c("control","accelerated snowmelt"), values = c("control"="darkgreen","manip"="blue"))
Again we plotted the pollinator activity period.
ggplot(pot,aes(x=as.Date.factor(pot$midpoint), y=pot$total.seeds, group= pot$treat)) +
labs(y="Total Developed Seeds",x="Date", title="Pollinator Activity Period and Total Developed Seeds in Potentilla") +
theme_classic() +
geom_point(size=3, alpha = 0.8, aes(color=treat)) +
geom_smooth(aes(color=treat)) +
scale_color_brewer(palette="Dark2") +
labs(color="Open or Hand-pollinated") +
scale_color_manual(labels= c("open","open-hand"), values = c("open"="darkred","open-hand"="darkorange")) +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
annotate("rect",xmin=as.Date.character("2019-07-05"), xmax=as.Date.character("2019-08-24"),ymin=0,ymax=Inf,alpha=0.1) +
scale_x_date(date_breaks = "3 days", date_labels = "%m/%d")
